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Abstract 

We present a detailed theoretical study of transverse exciton-polariton patterns in semicon- 
ductor quantum-well microcavities. These patterns are initiated by directional instabilities 
(driven mainly by polariton-polariton scattering) in the uniform pump-generated polariton 
field and are measured as optical patterns in a transverse plane in the far field. Based on a mi- 
croscopic many-particle theory, we investigate the spatio-temporal dynamics of the formation, 
selection, and optical control of these patterns. An emphasis is placed on a previously pro- 
posed low-intensity, all-optical switching scheme designed to exploit these instability-driven 
patterns. Simulations and detailed analyses of simplified and more physically transparent 
models are used. Two aspects of the problem are studied in detail. First, we study the de- 
pendencies of the stability behaviors of various patterns, as well as transition time scales, on 
parameters relevant to the switching action. These parameters are the degree of built-in az- 
imuthal anisotropy in the system and the switching (control) beam intensity. It is found that 
if the parameters are varied incrementally, the pattern system undergoes abrupt transitions 
at threshold parameter values which are accompanied by multiple-stability and hysteresis be- 
haviors. Moreover, during a real-time switching action, the transient dynamics of the system, 
in particular the transition time scale, may depend significantly on the proximity of unstable 
patterns. The second aspect is a classification and detailed analysis of the polariton scattering 
processes contributing to the pattern dynamics, giving us an understanding of the selection 
and control of patterns as results of these processes' intricate interplay. The crucial role played 
by the (relative) phases of the polariton amplitudes in determining the gains and/or losses of 
polariton densities in various momentum modes is highlighted. As a result of this analysis, an 
interpretation of the actions of the various processes in terms of concepts commonly used in 
classical pattern-forming systems is given. 
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I. INTRODUCTION 



Laser beams propagating through a nonlinear medium can under certain conditions 
undergo directional instabilities. These instabilities are driven by phase-conjugate wave- 
mixing processes,^® leading to spontaneous generation and build-up of intensity in 
modes with propagation directions other than those of the input beams. The system's 
transverse translational symmetry is thus spontaneously broken, resulting in periodic 
or quasiperiodic modulational patterns, e.g. stripes and hexagons, in the intensity pro- 
file. These spontaneously formed optical patterns have been observed in experiments 
with lasers,^^ and such passive optical systems as atomic gases J^^, liquid crystal 
light valves,^, photorefractive crystals,^ and, most recently, semiconductor quantum- 
well microcavitiesP^ (For general reviews, see e.g. Refs. [T51 - EU1 ) In a broader context, 
these developments are instances of spontaneous pattern formation found in a wide 
range of nonequilibrium physical systems such as classical fluids and chemical reaction 
systemsp3Z2 

The occurrence of directional instability patterns in optics has also indicated the 
potential of using them for low-intensity all-optical switching. It was demonstrated in 
Ref . HDthat the orientation of transverse patterns in a rubidium vapor can be all-optically 
and reversibly switched by the application of a weak control beam. Subsequently, it was 
proposed that a similar switching scheme can also be implemented in a semiconductor 
quantum-well microcavityP Another analysis also shows the possibility to optically con- 
trol transverse patterns in planar semiconductor structures.^ A review of these activities 
was given in Ref. [231 The switching scheme devised in Ref. [H] has recently been experi- 
mentally demonstrated in a quantum- well double-cavity system.^ In this context, it is 
worthwhile to note the different physical origins of pattern formation in these systems. 
In an optically excited atomic vaporp^ nonlinearities stem from partial saturation of the 
atomic resonances, and three-dimensional phase-matching conditions have to be satisfied 
with counterpropagating pump beams for efficient off-axis four-wave mixing of light to 
take place. In contrast, in a planar quantum-well based microcavity, nonlinearities arise 
from the parametric scattering of exciton-polaritons^^ (polariton four- wave mixing), 
which is predominantly driven by the Coulomb interaction between the polaritons' ex- 
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citonic constituents] 27 " 30 ^ with additional smaller contributions from phase-space filling. 
Two-dimensional phase-matching, in the cavity's plane, has to be fulfilled to render the 
scattering processes efficient. 

While setting the conceptual foundation of the directional switching action in semi- 
conductor microcavities in Ref. El and implementing the switching in Ref. [TT] were im- 
portant steps forward, further work is required to fully understand the underlying com- 
plicated dynamics. In this paper, we study this problem theoretically as an instance of 
competitions among various transverse patterns formed by the nonlinear interactions of 
exciton polaritons. The switching action of the control beam is viewed as destabiliz- 
ing the initial pattern while stabilizing another. Two aspects of this investigation are 
reported in this paper. First, we study the dependencies of the stability behaviors of var- 
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FIG. 1: (Color online.) (a) Sketch of the semiconductor quantum-well microcavity. Two 
distributed Bragg reflectors (DBR) form a photonic cavity with an embedded semiconductor 
quantum well (QW). The pump is in normal incidence to the quantum- well plane. When 
transverse instabilities are present, off-axis patterns can be observed in the far field. The 
orientations of these patterns can be rotated by applying a weak control beam, (b) Cavity 
polariton energy vs in-plane momentum. The uncoupled exciton and cavity photon modes are 
shown as dotted and dashed lines, respectively. The solid lines represent the upper (UPB) 
and the lower (LPB) polariton branches. The displayed dispersion relations do not include 
the effects of (nonlinear) exciton interactions. The pump frequency and scattering of two 
pump-induced polaritons to off-axis modes are also indicated. 
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ious patterns, as well as transition time scales, on parameters relevant to the switching 
action. These parameters are the degree of built-in azimuthal anisotropy in the system 
and the control beam intensity. Second, we analyze in detail the actions of physical 
processes contributing to the dynamics and form an understanding of the selection and 
control of patterns as results of these processes' intricate interplay. 

The nonlinear optical and many-boson aspects of polaritons in semiconductor mi- 
crocavities are an active area of current research.^ Besides modulational intensity pat- 
terns, more complicated spatial instability structures such as solitons^l and vortices®' 
are also being investigated in these systems. Optical parametric oscillation in various 
pump/signal configurations have been demonstrated in quantum well microcavitiesP^^ 
This process shares a common underlying mechanism - parametric polariton scattering 
- with the directional instabilities studied here. 

The quantum well microcavity configuration used in this paper is sketched in Fig.[T](a). 
A photonic cavity is formed by a pair of distributed Bragg reflectors (DBRs), and a semi- 
conductor quantum well (QW) is placed inside the cavity. One cavity photon mode is 
tuned to coincide with the lowest Is heavy hole excitonic resonance, leading to strong 



1 1 




4 4 



FIG. 2: (Color online.) Schematic representation of the switching process. The sites labeled 
1 to 6 represent the projected positions on a transverse plane in the far field of six off-axis 
momentum modes, forming a hexagon. On the left is shown an optical pattern with the mode 
pair 1 and 4 being spontaneously generated by phase-conjugate instabilities and favored by a 
built-in anisotropy of the system. When a control beam is directed at mode 2, the pattern is 
rotated to directions 2 and 5. When the control is then turned off, the pattern reverts to its 
original orientation. 
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coupling between the two and the formation of an upper (UPB) and a lower (LPB) po- 
lariton branch. The dispersion relations (energy vs in-plane momentum) of these linear 
optical excitations for the parameters used in this paper are plotted in Fig. [TJb). A 
pump beam comes in at normal incidence to the cavity's plane and is spectrally tuned 
above the bottom of the lower polariton branch but well below the bare exciton reso- 
nance. It excites (virtual) polaritons with zero in-plane momenta, pairs of which can 
scatter near-resonantly into modes with finite, opposite transverse (in-plane) momenta 
on the LPB, as indicated in Fig. [I](b). The scattered polaritons have a certain prob- 
ability to emerge from the microcavity as photons. Under favorable conditions, when 
the pump intensity is above a threshold, the emerging photons in the off-axis directions 
become intense and coherent beams, characteristic of parametric, directional instabili- 
ties of the pump beam. The instability reduces the system's symmetry from transverse 
translational to azimuthal. This azimuthal symmetry may be further broken by the 
nonlinear interactions among the off- axis polariton field amplitudes, thereby favoring 
patterns with density in a finite number of clusters of k modes (i.e. modes in transverse 
momentum space), e.g. a hexagon or its subsets. Schematically illustrated in Fig. [2j 
the low-intensity switching scheme proposed in Ref. |H] involves transitions between two 
two-spot patterns residing on the same hexagon. (Note the role of the built-in azimuthal 
anisotropy in the preparation of the initial pattern.) Instead of coherent polariton scat- 
tering, the instabilities may equivalently be interpreted as arising from phase-conjugate 
four wave mixing of the coherent polariton field. 

The conditions for the formation of hexagons and two-spot patterns are quite well 
understood in other disciplines of physics.^ However, the additional degrees of complex- 
ity, anisotropy and external real-time control, make pattern formation in our system, 
and in the system in Ref. HH a much more involved issue than found in commonly 
studied systems. Moreover, unlike macroscopic pattern-forming systems, where the field 
is real-valued, our polariton field is quantum mechanical and therefore complex-valued, 
and the mutual activating and inhibiting actions of the mode densities on each other are 
entwined with their relative phases. This makes building physical intuition about these 
actions and hence how patterns are selected a non-trivial exercise. 
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Our theoretical approach is based on a microscopic theory of excitonic optical re- 
sponse of the semiconductor quantum wel l 27 ' 28 1 coupled to a classical treatment of the 
light field. Only pairwise interactions among excitons, evaluated in the Hartree-Fock 
(HF) approximation, and Pauli blocking modifications of the light-exciton coupling are 
included in the nonlinear exciton dynamics. The theory is discussed in more detail in 
Section [IT] below. To facilitate conceptual understanding, we treat only the simple case 
of a circularly polarized polariton field in a single quantum well cavity [Fig. [I]. Ow- 
ing to the spin dependence of the exciton-exciton interaction,^ parametric scattering 
for linearly polarized polaritons is different from that in the circularly polarized caseP^ 
There are indication d 17 * 35 * 37 -* that transverse instabilities are more favored in the linearly 
polarized case. Moreover, using a double cavity also makes pattern formation easier.^ 
We believe the conclusions drawn in this paper about physical mechanisms also apply 
to these more complicated cases. 

We solve the equations of our theory in various model state spaces. At the most 
complete level, we solve the real(r)-space version of our equations directly on a grid. 
While such simulations are realistic, they are computationally demanding and it is not 
easy to disentangle the roles played by various dynamical processes in generating the 
results. For this latter purpose, we also evaluate the theory within reduced model 
state spaces consisting of small numbers of k modes, the relative simplicity of which 
makes them much more amenable to analysis. In our parameter-dependency study and 
our analysis of the contributing dynamical processes, we make heavy use of a single- 
hexagon model which has a state space consisting of one pump mode (k = 0) and six 
off-axis k modes arranged in a regular hexagon. Larger models are also used to analyze 
competitions among hexagons and more extended patterns. 

This paragraph and the next summarize our findings. For typical GaAs parameters, 
even if the system has complete azimuthal isotropy, a ring pattern in momentum space is 
destabilized in favor of a single hexagon with arbitrary orientation. A small anisotropy is 
sufficient to destabilize the hexagon in favor of a two-spot pattern. As one increases the 
anisotropy continuously from zero, the transition from hexagon to two-spot is abrupt, 
and below this transition, a range of anisotropy values exist where both patterns are 
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stable. In the switching scheme, a two-spot pattern stabilized by anisotropy is set up as 
the initial signal, and the control beam is directed at another mode pair (see Fig. [2]). For 
our chosen material (GaAs) and pump beam parameters, the pattern can be switched 
by a control beam of intensity 1/6500 times the intensity of the two-spot signal. In a 
detailed parameter variation study, when the control intensity is increased from zero, 
the stable state transits from the two-spot pattern favored by the anisotropy, through 
an asymmetric hexagon (i.e. a hexagon with asymmetric distribution of intensity), to 
the two-spot pattern favored by the control. The transitions are again abrupt, with a 
bistable region at each transition. The switching time scale is normally of the order of 1 
ns but diverges as the control intensity approaches from above the higher critical value 
for transition. 

The analysis of the mechanisms underlying these phenomena has yielded the following 
picture. The dynamics of the polariton amplitude in a particular (off-axis k) mode is 
driven by pairwise scatterings for which the mode in question is one of the outgoing 
modes. We classify the processes as linear, quadratic, or cubic according to how many 



of the other three modes involved in the scattering are off-axis (Fig. 14). The linear gain 
processes provide the basic instability mechanism to the pump-excited uniform polariton 
field. For a pump intensity slightly above threhold, the set of k modes experiencing 
linear instability growth form a ring of finite thickness around the origin in (in-plane) 
momentum space. As the polariton amplitudes grow in these modes, the quadratic and 
cubic scattering processes, which represent interactions among the off-axis modes, come 
into play. They drive the competitions among phase-conjugate mode pairs (pairs with 
equal but opposite momenta (k and — k)) which in turn determines whether a particular 
pattern (ring, hexagon, two-spot etc) is stable. Because of momentum conservation, the 
quadratic scatterings take place only among off-axis modes that are on a regular hexagon 



(see Fig. 14 below), whereas the cubic processes are not subjected to this restriction. 
Therefore, competitions among modes residing on different hexagons are driven only by 
the cubic processes. It is important to note that the scatterings occur among coherent 
amplitudes: the relative phases of the involved polariton amplitudes play a crucial role 
in determining whether a particular process contributes to a gain or a loss in density for 
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the mode in question. Because of the complex feedback between the dynamics of the 
phase and the magnitude of the polariton field, it is difficult to cleanly disentangle the 
gain/loss effects of the different processes. Even so, from the simulation results one can 
still characterize their effects using concepts common in classical pattern forming systems 
as follows. For the typical GaAs parameters used here, the stability of the hexagonal 
pattern over two-spot patterns (in the isotropic case) indicates that the net effect of the 
quadratic processes is 'cross-activating', i.e. the presence of density in one conjugate 
mode pair favors the increase of density in other mode pairs. The main effect of the 
cubic processes is to bring about saturation of the off-axis polariton density. One can 
subdivide them into 'self-saturating' (the density in a conjugate mode pair inhibiting its 
own growth) and 'cross-saturating' (the density in one mode pair inhibiting the growth of 
density in other pairs) processes. The spontaneous breaking of the azimuthal symmetry 
of the system, favoring the radially localized hexagons over a ring pattern, indicates that 
the cross-saturating processes are stronger than the self-saturating ones. 

The above findings are based mostly on analyses of models in which the pump in- 
tensity is infinitely localized in momentum space (in the mode k = 0), or equivalently, 
uniformly distributed over the 2D real space of the cavity's plane, and the state space is 
truncated. Further investigations into the effects of a finite pump laser spot size and in- 
clusion of more modes in the state space will be carried out within the full 2D real-space 
simulation model. 

The paper is organized as follows: We first describe and discuss our theoretical ap- 



proach in Section [11} In Section III, we show a set of results from a 2D r-space simulation 



and discuss two reduced models in truncated sets of k modes. In Section [TV] we introduce 
and validate the single-hexagon model through comparison with the previous models. In 
Section |Vj we present the effects of variations of anisotropy and control beam intensity 
on pattern selection and switching behavior in the single- hexagon model. In Section 



VI, we classify the scattering/four wave mixing processes among the polariton modes, 



and the complex interplay of these processes are analyzed in Section VII| An analysis of 
the competitions among modes belonging to different hexagons is carried out in Section 



VIII Section IX contains some concluding comments. 
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II. THEORY OF NONLINEAR OPTICAL RESPONSE OF A QUANTUM 
WELL MICROCAVITY 



This section outlines the microscopic theory of the third-order nonlinear response of 
a quantum well microcavity, from which our working equations are derived. The theory 
can conceptually be divided into two parts: the nonlinear response of the quantum well 
to the presence of cavity photons and the coupling between the relevant cavity photon 
modes and the light field outside of the cavity. Since we work in the regime of strong 
linear coupling between the QW exciton and a cavity photon mode, it is beneficial to 
consider the polaritons as elementary excitations of the QW microcavity and the cavity's 
nonlinear response as being due to the interactions among the polaritons. For the theory 
of nonlinear QW optical response, a more detailed account can be found in e.g. Refs. 
27,28,38,39. An account of our present theory was given in Ref. |23] 

In Ref. |23l our treatment of the passage of light into and out of the cavity, when a 
given polaritonic excitation is present inside the cavity, was only very briefly explained. 
We start this section with a fuller account of this part of our theory. Numerically, the 
most accurate way to propagate light waves across the microcavity is, of course, to 
solve the Maxwell equations throughout the DBR+QW structure using e.g. a transfer 
matrix or a finite-difference time-domain method. 42 However, since we 

are concerned with obtaining an overall picture of parametric dependencies rather than 
simulating a specific experimental setup, we have chosen to use a simple model that 
approximates the wave propagation effects as a coupling between the radiation field out- 
side of the cavity and the chosen cavity mode. We will call the former the 'macroscopic 
field'. In this model, the whole microcavity is treated as a zero-thickness structure which 
confines both the cavity photon field and the exciton field in the longitudinal direction 
and lets them extend across the entire transverse plane. The effect of the cavity pho- 
ton field on the macroscopic light field is represented as a 'polarization density' in the 
wave equation governing the latter. Explicitly, we write the electric displacement of the 
macroscopic field (which includes the incident, reflected, and transmitted waves) as 

V(r, t) = e n 2 s £(r, t) - M c E cav (x, y, t)5(z) (1) 
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where S(r,t) is the macroscopic electric field, eo is the vacuum permittivity, n s is the 
refractive index of the substrate outside of the cavity, E cscv (x,y,t) is the cavity photon 
field, and t c is a coupling strength. For later convenience, the dimension of Ht c is chosen 
to be that of electric charge. In Eq. we have set up a coordinate system in which 
the z-axis is along the longitudinal direction and the cavity is positioned at z = 0. 

The incident waves (pump and control) come in from the left and are all '+ '-circularly 
polarized. The exciton spin states in the QW are quantized along the z-axis so that an 
obliquely incident (+)-polarized beam may also generate excitons with spin= —1 which, 
after interacting with other excitons, may generate a (— )-polarized component in the 
outgoing waves. However, the polar angles of the off-axis beams are usually quite small, 
and we, for simplicity, ignore this complication and take all outgoing fields to be (+)- 
polarized. The algebraic development within the model is detailed in Appendix |Aj with 
the following results. Under the assumption that each of the waves involved is a slowly 
varying envelope modulating a plane wave with a common carrier frequency u p , the 
electric field in the half-space z < can be written in the form 

S(r,t) = J2e l{kxX+kyy) e + [E Kmc (t + k z (k,uj p )z/io p ) - E^ mR (t - k z (k } oo p )z/uj p )} ,z<0 

k 

(2) 

where k = (k x ,k y ), k = |k|, and k z (k,u p ) = + uj 2 n 2 / c 2 — k 2 , c is the light speed in 
vacuum, and e + is the unit helicity vector for the (+)-polarized field. For each transverse 
momentum k, -Ek,inc is the incident field and i?k,refl the reflected field. In the z > half 
plane, the field can be written as 

S(r, t) = J2 e i[kxX+kvy) e+E k , tIWS (t - k z (k, u p )z/u p ), z>0 (3) 

k 

where i?k,trans is the transmitted field. The cavity field, which by our argument above is 
also (+)-polarized, is also decomposed into its spatial Fourier modes: 

E cav (x, y,t) = J2 e i{k * x+k * y) e+E k (t) (4) 

k 

We stress that E k (t) here denotes a momentum-space component of the cavity field and 
is not to be confused with the macroscopic field Sir, t). The k sums in Eqs. ([2])-((4]) range 
in principle over the whole (transverse) momentum space, but they are restricted to a 
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finite set of k modes in the reduced-mode model calculations discussed in this paper. 
For each k mode, the transmitted (-Ek.trans), and reflected (-Ek.refl) light fields are solved 
in terms of the incident field (.E^mc) and the cavity field (E^) in Appendix [Aj The result 
is: 

-^k.trans = -^k.inc — -^k,refl (5) 

ht c dEk 

E ^ cR = ~2^^T (6) 

We next postulate the equation of motion of the cavity field under the actions of the 
macroscopic light field and the exciton field. In writing down the equation, we make the 
reasonable assumption that each transverse k mode of the cavity field, E^, is a simple 
oscillator driven linearly by the other two fields. The proportionality coefficients in the 
coupling terms are chosen in such a way that \E^\ 2 can be interpreted as mode photon 
density (see Appendix IaI). The equation is: 



where cj£ is the uncoupled mode dispersion relation for a cylindrically symmetric cav- 



ity: w£ = a/wq 2 + c 2 k 2 /n 2 , Uq being the cavity mode energy at k = 0, and n& being 
an effective refractive index of the cavity's medium which characterizes the transverse 
propagation speed of the cavity mode, pk denotes the exciton field, which will be identi- 
fied below with the Is heavy-hole exciton contribution to the interband polarization in 
the QW, so that |pk| 2 is the exciton density (number of excitons per unit area) in the k 
modeP^ The exciton-photon coupling strength is an input parameter in Eq. ([7| but 
will be given a microscopic meaning in terms of electron-hole dynamics below. E^ inc is 
the (limiting value of the) macroscopic electric field at the cavity's position. Within our 
model, we have 

^Lc(*) = -^k.transW (8) 

The coupling coefficient t c was introduced in Eq. ([I]) above. If we set the incident field 
-Ek.inc and the exciton field pk to zero, Eqs. ([5])-(j8]) would reduce to a homogeneous 
equation for E^ from which we can extract a radiative decay rate (in energy units) 
T = ujph 2 tl/(2eocn s ) for the cavity. We reiterate that Eq. Q is written down as a 
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postulate based on physically reasonable assumptions. This treatment is sufficient for 
our purposes here. Alternatively, it could be derived as a classical approximation in a 
theory that starts with the cavity photon field being quantized. 

We next briefly describe the microscopic theory of the excitonic optical response of 
quantum wells which gives the equation of motion of p k l 27 l 28 l 3 S B9] starts with a many- 
particle Hamiltonian operating on the Fock space of electrons and holes. Referring 
to Fig. [T] again, for the quantum well, we consider a GaAs-type semiconductor band 
structure around the fundamental band gap with an exciton size of about 10~ 6 cm. For 
the present study, the single-particle basis includes only one parabolic conduction band 
and one parabolic valence band, i.e. the heavy hole band. The Hamiltonian contains 
the charge carriers' kinetic energy and pairwise Coulomb interactions among them, and 
a coupling to the cavity photon field that either creates or annihilates an electron-hole 
pair in the rotating wave approximation. Thus the only material parameters in the 
theory are the electron mass, the hole mass, the background dielectric constant, and the 
interband dipole coupling. With this Hamiltonian, equations of motion of observables 
(certain expectation values of products of charge carrier creation/annihilation operators) 
are derived. The dynamical variable relevant to us, the exciton field Pk{t), is defined 
in this microscopic theory as the interband polarization restricted to the heavy-hole Is 
exciton subspaceP^ 

= 72 S^ k ' + WK-ic(*Kic+k(t)) (9) 

k' 

where a e ,q(t) is the Heisenberg-picture annihilation operator of a conduction band elec- 
tron with momentum q at time t and a/ l . q (t) is the corresponding operator for a hole 
state, and (• • •) represents the expectation value in the initial state of the electron-hole 
system. /3 = m,h/(m e + m/J is the ratio of the hole mass to the total exciton mass, 
C 2 is the area of the normalization box, and 0(q) = v^r^o/fl + (ao<?/2) 2 ] 3//2 is the 
two-dimensional internal Is exciton wavefunction in the space of relative electron-hole 
momentum, with ao = h 2 €b/ (e 2 m r ), e& and m r being the (3D) exciton Bohr radius, the 
background dielectric constant and the reduced mass of the electron-hole pair respec- 
tively. The restriction to the Is state is valid because we limit the pump intensity so 
that the maximum generated polar iton density is of the order of 10 10 cm~ 2 . At this 
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density and the chosen pump frequency, which is tuned to a frequency within the lower 
polariton branch (Fig. [lib)), the excitations in the quantum well are expected to stay 
predominantly in the Is heavy hole exciton state. For a more detailed analysis on the 
validity of the Is approximation, see Ref. |28j Since we restrict our considerations to 
+-polarized excitons, the electron and the hole in Eq. ^ have spins —1/2 and 3/2 
respectively. 

Starting from the electron-hole vacuum as the initial state, the equation of motion for 
the Pk(t) is derived up to third order in the light field amplitud e 1 38 1 39 1 43 1 on the Hartree- 
Fock- level. In the derivation, the single-pair electron- hole basis is restricted to the Is 
state as mentioned. In this limit, the equation for a system of co-circularly-polarized 
excitons in the coherent third-order regime isj 8 * 27 * 28 ! 

^'dT = ^ _ ^ x ">P k ~ ^ k ^ k + ?2 ( 2 ^k''£>q£>k'-Ek'' + VHFPqPk'Pk") ^q,k'+k»-k 

qk'k" 

(10) 

Here e£ is the free exciton energy and 7^ is a phenomenological dephasing rate which 
represents all non-radiative losses, including the transfer of coherent excitons to the 
incoherent population. The photon-exciton coupling Q^, which is an input parameter in 
the cavity photon equation, Eq. ([7]), is given in the present microscopic theory by (in 
S.I. units) 



^k = e(r) eh 0*(O)/ 7 (2 QW )y M,/e e 6 , (11) 

where (r) e h is the interband dipole matrix element, 0*(O) = 2-y/2/(ao\/7r) is the two- 
dimensional real space Is exciton wavefunction (Fourier transform of 0(q)) at zero 
electron-hole separation, f^z) is the resonant cavity photon wavefunction along the coor- 
dinate axis (z-axis) normal to the cavity's plane, normalized to unity: J™ dz\f 1 (z) | 2 = 1 
and zqw is the position of the quantum well on the z-axis. 

Of the two third-order terms, the term with the pp*E structure stems from Pauli 
blocking, or phase-space filling (PSF), among the fermionic components of the exciton 
that are created by photon absorption and those of excitons already in the system. The 
strength is given by A = 2na^jl . The other term is due to mean-field, or Hartree- 
Fock (HF), Coulomb interaction among the excitons, with strength Vhf = I^E^a^iX — 
3157r 2 /4096) ~ 1.52E b al, where E b is the (2D) exciton Rydberg energy. Both A and 
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Vhf actually depend on the momenta of the involved excitons and photons, but since 
the dependencies are rather weak and the relevant momenta are not excessively large, 
we approximate them by their values at zero momentum. Each exciton interaction or 
PSF term can be visualized as either a four-wave mixing process with wave-vector and 
frequency matching or a polariton-polariton scattering process with momentum and 
energy conservation. The full set of third-order terms also includes a time-retarded 
Coulomb correlation between two excitons, which is neglected here. Because of its time- 



nonlocal structure, including the Coulomb correlation term in Eq. (10) would increase 
the complexity in solving the equation substantially. It was previously showrpZES that it 
is much less important than the PSF and HF terms when the two interacting excitons are 
co-circularly polarized and have energies far below the exciton resonance. So dropping 
the retarded- Coulomb term in our case is expected to be a good approximation. A 
detailed study of its effects on polariton pattern formation will be a subject for future 
work. 

Moreover the dynamical effects of the dephased excitons are not considered. A more 
comprehensive theory would include a description of the transfer by dephasing and the 
subsequent interactions between the coherent and incoherent populations. Equations 
of motion serving this purpose could be derived within the microscopic frameworks 



that produce Eq. (10). While recent experiments indicate that the incoherent exciton 
population remain smalP^ in the density and frequency regions that we consider, its 
cumulative effects should be studied more carefully. Similar considerations apply to 
three-exciton and larger-cluster correlations which would be accounted for in an extended 
theory beyond the regime. 

To recapitulate, Eqs. ([5]), ([6]), dTl), and (10) are the working equations in the 
present paper. With the incident macroscopic light field Ek,i n c(t) as input, the response of 
the quantum well microcavity, E^t) and Pk(t), as well as the resulting output light fields, 
-E-k.transW an d -E'k,refl(^), are calculated by solving this set of equations simultaneously. 
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III. SIMULATIONS AND REDUCED MODELS 



In this section, we show some representative numerical results on transverse pattern 
competition and control in a quantum well microcavity by solving, at various levels of 
approximation, Eqs. ([5])-(|8]), and (10) laid down in the previous section. The numerical 
challenge stems mainly from the nonlinear (HF and PSF) terms in the exciton field 



equation, Eq. (10), which couple the k modes together. Since the planar normalization 



box (in real space) is taken to infinity, the summation sign in Eq. ( 10 ) actually denotes 
k space integrations. A comprehensive simulation of the dynamics of the quantum-well 
microcavity then involves solving the equations on a two-dimensional momentum space 
grid, or, alternatively, solving the Fourier transforms of the equations on a real space grid 
in the plane of the microcavity. While such simulations are realistic, they are computa- 
tionally demanding, and their results are not easy to analyze. Since, when the system 
is not far above the instability threshold, the stationary transverse patterns are usually 
distributed over localized regions in k space, approximations in which the equations 
are solved within judiciously chosen, vastly reduced sets of k modes can profitably be 
used. In this paper, we take advantage of this fact and put our emphasis on dissecting 
the results from several such reduced models in order to advance towards our aim of 
understanding the physical mechanisms of pattern competition in our system. 

In the following, after some preliminary remarks on the linear instability of the spa- 
tially planar homogeneous (with only the k = component) state, we will present a 
set of representative results of our full two-dimensional simulations. We will then intro- 
duce two reduced models - the 'multi-|k| model' and the 'ring model' - in which the 
equations of motion are solved over selected sets of k modes. These two models essen- 
tially reproduce the representative results of the full 2D simulations. We will extend 
our investigations using these models. Based on these investigations, we will argue at 
the end of this section that the analysis of an even simpler model, the state space of 
which consists of the k = mode and a hexagonal set of six off-axis modes, could be 
expected to provide the critical insights to understand the patterns' behaviors in the 
more comprehensive models. This 'single- hexagon' model will be analyzed in detail in 
the following sections. 
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For the calculations, material parameters appropriate for GaAs are used: Ey, ~ 
13 meV, a w 170 A, e x = 1.497 eV, ^ x = 0.4 meV, Q k = Q = 8meV, n b = 3.6, 
and eb = n\. With the given values of Ej, and ao, Vhf ~ 0.57 x 10 _10 meVcm 2 , and 2AQ 
is about three times weaker. For the cavity mode, we choose Fuvq = Eq, and t c is chosen 
such that the radiative decay rate is T f» 1.5 meV for Hu p around 1.5eV. The incident 
pump is tuned to 5 meV below the bare exciton resonance, i.e. hw p = 1.492eV, and its 
steady state flux intensity, which we denote by Jp mp , is about 25.2 kWcm -2 . 

at the lossless limit, i.e. with Vhf, A, j x , and 
t c set to zero, give the linear polariton spectrum and the polariton eigenvectors in the 
photon-exciton basis. With the other parameters set as above, the linear polariton 
spectrum is plotted in Fig. [TJ Returning to the full nonlinear equations, one can easily 
verify that when the input field is incident normally on the cavity (pump only), there 
exist solution(s) that preserve the planar translational symmetry of the input beam (pk 
and E\f_ vanish for all k ^ 0). At steady state, po and E$ oscillate at the pump frequency 
Up. The exciton density \po\ 2 calculated for the above parameter values is plotted as 
a function of the incident field intensity I{ nc = |-Eo,inc| 2 i n Fig. |3j A linear stability 
analysis with off-axis (k ^ 0) perturbations on the uniform solution is performed, and 
the range of instability against this class of perturbations is also shown in Fig. [3j The 
stability matrix is block-diagonal in k space with the fields pk, -Ek,P-k, and ELk for 

6n 
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FIG. 3: (Color online.) Steady state pump-generated exciton density in the quantum well as a 
function of the pump beam intensity. For densities inside the red, dashed segment, the uniform 
polariton field is unstable against finite |k| perturbations. The cross marks the pump intensity 
used in the simulations here. 



The linearized Eq. <\7h and Eq. (10 
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each k coupled together. For each value of |po| within the off-axis instability range, 
the uniform state is unstable against small random perturbations in a set of k modes, 
which form a ring of finite radial thickness in k space. The mean radius of this ring is 
approximately determined by phase matching that fulfills the resonance in the polariton 
scattering, as illustrated in Fig. ^p. The dispersion relations in Fig. [T|d are calculated 



without the exciton interactions in Eq. (10), which shift the exciton energy, and hence 



the (lower) polariton energy, upwards. As a result, the actual resonant |k| value is 
smaller than that shown in Fig. [TJd. The point (|-£/o,inc|> \po\) corresponding to the 
incident field strength used in our calculations is marked in Fig. |3j For this value, the 
instability-triggering k modes lie in the range 1.25 x 10~ 4 cm _1 < |k| < 1.40 x 10 _4 cm _1 . 
More detailed discussions of the linear stability analysis of the uniform state may be 
found in Refs. I3TH41 In our calculations, starting from the uniform state under steady 
pump irradiation, the polariton field (the lower polariton branch) for each mode inside 
the 'ring of instability' initially grows in time exponentially from random fluctuations. 
Subsequently, the interactions among these off-axis polariton modes, and the k = 
mode, govern the long-time competition dynamics among them. 

A. 2D real-space simulations 

The full 2D simulations are performed in real space by solving the Fourier transforms 
of Eqs. ([7]) and ( 10 ). We show a representative set of results here. A more comprehensive 



discussion of these simulations will be presented in a future publication. The configu- 
ration space is a square (spatial) grid with a step size of Ax = Ay = 0.45/zm in a box 
of length 90/im on one side. The initial state of the system is the polariton 'vacuum', 
i.e., the ground state of the semiconductor microcavity. The pump beam is switched 
on at t — 0, its intensity becoming steady after 8 ns. Shown in Fig. |4|a), the steady- 
state pump profile is broad and flat so that its corresponding (transverse) momentum 
distribution is narrowly peaked at k = 0. As time progresses, coherent polariton am- 
plitudes in several off-axis k modes form "spontaneously", causing density patterning 
inside the cavity. The calculated real-space exciton density distribution inside the QW 
at t = 18ns is shown in Fig. [4](b) . The corresponding momentum distribution is plotted 
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FIG. 4: (Color online.) Two-dimensional real-space simulations of hexagonal pattern for- 
mation, (a) The pump profile in 2-dimensional real space, (b) The exciton density profile in 
the QW with modulations stemming from transverse instabilities, (c) Hexagonal pattern in 
transverse momentum space formed by the exciton field in (b). The density in the center (in 
the pump mode) has been masked for clarity. See text for details of the simulations. 

in Figure ilc), which shows a dominant hexagonal pattern situated within the ring of 
linear instability. The azimuthal orientation of the hexagon is random. 

B. Reduced models 

The linear off-axis instability spontaneously breaks the (transverse) translational sym- 
metry of the system's setup but preserves the latter's azimuthal symmetry. In the full 
2D simulations, when the off-axis polariton fields are being built up, their mutual inter- 
actions also break the azimuthal symmetry, resulting in a stable hexagonal structure. 
As explained above, the simulation results indicate that the collection of modes taking 
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part in the pattern dynamics are relatively localized in k space. We take advantage 
of this fact and use reduced models in the following to investigate the competitions 
among different patterns, their control, and the physical mechanisms underlying these 
competitions. 

In the reduced models, we restrict our state space to a selected finite set of k modes. 
We introduce two such models in this subsection: (i) The 'multi-|k| model'. For the 
state space in this model, we choose six directions corresponding to the vertices of a 
regular hexagon centered at k = 0. The hexagon's orientation is arbitrary. N evenly 
spaced points along a radial segment in each direction are included in the state space, 
as illustrated in Fig. p](a). Explicitly, we write the k points as: 

K,i = (ko + h5k)ei , h = 1, • • • , N , i = 1, • • • , 6 (12) 

where e iy i = 1, ■ ■ ■ ,6 are unit vectors in the vertex directions, and the minimum radius 
k and the grid size 5k are run-time parameters in the numerical calculations, (ii) The 
'ring model'. In this model, we choose the |k| modes to lie on a circle centered at the 



origin, regularly spaced in the angle coordinate, as illustrated in Fig. |5[b). We also 
organize the points into iV hexagons, writing them as 



k r e h)i , h = 1, • • • , N , i = 1, • • • , 6 (13) 



where e^j is a unit vector in an arbitrarily selected direction, eh+i,i is obtained by 
rotating e^j counterclockwise through the angle and the radius of the circle, k r , is 
a run-time parameter. 

By considering the momentum conservation constraint in the nonlinear terms in 



Eq. (10), one can see that, for both models, Eq. (10) reduces to equations of the follow- 
ing form (we have simplified the notation pk h ; to p hi and po to po). The k = (on axis) 
equation reads: 
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FIG. 5: (Color online.) Sketches of the state spaces, made up of modes in transverse momentum 
space, of the three reduced models used in this paper, (a) Multi-|k| model. Beside the origin, 
the modes lie on six radial segments (thick dashed lines) arranged hexagonally. (b) Ring model. 
The off-axis part of the state space is a ring of linearly unstable k modes, (c) Single-hexagon 
model. The state space is the origin plus a six-mode subset of the ring model state space that 
form a regular hexagon. 
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The off-axis equations, for h = 1, • • • , N and i = 1, 



, 6, read: 



ihp h 



-n 



+ ^HF 



ij x + 2V RF 



2A \ Po 



( hi 2 - \ \Ph,i\ 2 + ^2\Ph',i'\ 2 ) + 2A ~n [p*o E o + ^2Ph>,i> E h>,i> 

\ " h',i> J \ h',V 

- \Ph,i\ 2 + J2\Ph'A 2 1 



PoE + 22,Ph',i' E h',i'^ ~ {Ph,i E h,i+3 + Ph,i-3 E h,i) 
h',i' 

lPh',i-2 + Ph',i+2Ph',i-l) + Ph',iPh',i-3 \ 

h' h'^h J J 



+2P0 



AM (ph+2 E h,i+l+Ph,i-2 E h,i-l) + V HF {P*h,i+2Ph,i+l+Ph,i-2Ph,i-l) 



+2AQ (pl i+2 Ph,i+i + Ph,i-2Ph,i-i) E o 

When the off-axis direction number subscript i is outside the range 1 — 6, it represents 
the direction either i — 6 for i greater than 6 or direction i + 6 for i smaller than 1. 
For simplicity we have neglected the k dependence of and replaced it by fi — fio- 
The cavity field equations, Eqs. (J5|)- (J8|) , remain the same for each k hii , with analogous 
notations being used for the cavity and macroscopic fields: Ehi etc. 



On the right hand side of Eq. (14), the terms represent, in the order of their 
appearance, the HF shift in the exciton energy, photon scattering off an exciton density 
grating, Pauli blocking of the photon-exciton coupling, and various density transferring 



scattering processes. The terms in Eq. (15) can also be interpreted the same way. We 
will discuss the density transferring processes in more details below. 



1. The multi-\k\ model: results 

Within this model we explore the competitions among patterns that exist in a hexag- 
onal geometry and the radial distributions of these patterns in k space. The state space 
consists of iV = 256 |k| points in the direction of each of the six hexagonal vertices. 
We use a grid size of 22.5 cm -1 along each radial segment so that the state space is 
spread over the range from 1.10675 x 10 4 cm" 1 to 1.6805 x 10 4 cm -1 in each segment. 
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In this state space, Eqs. Q-Q, Q, and g5|) are solved in the rotating frame of 
the incident field, i.e. with e~ tuJpt factored out of all the fields in the equations. The 
time-stepping is performed with a 4-th order Runge-Kutta method with a time step of 
20 fs. This time step is sufficient for simulating the system dynamics, which has a typical 
time scale of the order of 1 ps, with good accuracy. Random fluctuating light sources 
(1.26 x 10 -6 Wcm -2 ) are added in the off-axis modes as seeds for the initial instability 
growth. 

Solving this model with the same material and beam parameters as in the full-2D 
simulations yields similar results, i.e., starting with the 2D uniform state, the system 
spontaneously generates off-axis coherent fields which subsequently stabilize as a hexag- 
onal pattern. Along each radial segment, the off-axis field initially grows in modes 
over a range of |k| values. When the pattern is stabilizing, however, one mode in each 
hexagonal direction wins, the intensity in the other modes in the segment decaying to 
zero. 

In the experiment using atomic vapors in Ref. [HI a stable two-spot pattern was 
obtained instead of a hexagon. It was theorized that an imperfect azimuthal symme- 
try, which is unavoidable in practice, might have favored the two-spot pattern. In the 
theoretical calculations on microcavities in Ref. El it was found that a small anisotropy 
introduced in one orientation is sufficient to favor the two-spot over the more symmetric 
distribution over all six hexagonal directions, and the subsequent controlled switching 
involves another two-spot pattern. For a recent detailed discussion of anisotropy in a 
real system, see Ref. US 

We have carried out calculations for this slightly anisotropic situation. The anisotropy 
is imposed, as in Ref. El by lowering i>j c k by 0.12 meV in directions 1 and 4 at every |k|. 
The values of other parameters remain as set above. A set of numerical results that 
include reversible switching by a control beam are shown in Fig. [6} The top three panels 
show the light intensities reflected from the microcavity in the six off-axis directions, 
and the bottom panel shows the control beam, as functions of time. In each direction, 
the intensity has been summed over all the modes on the radial segment in Fig. |5](a). 
The simulation starts at time t = with all fields equal to zero except for some low-level 
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fluctuations. The external pump field is turned on in the k = mode and reaches 
its steady state at about 30 ps with \p \ 2 = 2.6 x 10 10 cm -2 . The pump steady state 
intensity is just above the instability threshold, leading to the build-up of field intensity 
in the off-axis modes. The build-up occurs in all six directions initially, but after a while, 
the asymmetry (anisotropy) introduced in the cavity mode's dispersion relation gives a 
decisive advantage to directions 1 and 4, where the field intensity continues to grow to 
a steady state at around t ~ 2.7 ns, while the intensities in the other four directions 
fall back to fluctuation level. A control beam of intensity = 9.1 x 10~ 3 Wcm" 2 = 
3.6 x 10~ 7 ip Ump is introduced in direction 2 (see below for the radial |k| value of the 
control) at t = 3 ns, and it switches the off-axis intensity to directions 2 and 5, reaching 
a steady state level after 1.6 ns. In the switching process, the instability-generated field 
intensity is about 6500 times stronger than the control intensity. The control beam 
is then turned off at t = 5.5 ns, whereupon the signal reverts to directions 1 and 4, 
reaching steady state at t — 7.0 ns. This reversible switching process can be repeated 
indefinitely, as shown. We recall that the energy-scale parameters, i.e. VnFbo| 2 , 1x etc, 
in the polariton equations Eqs. ([7]), ( |l4| ), and (15) are of order of meV, giving a natural 
dynamical time scale of ps. However, the (linear) gain and loss rates in the off-axis 
modes offset each other by roughly two orders of magnitude (the pump intensity level 
is set here slightly above threshold) so that the changes in polariton density happen on 
a longer time scale between several tens of ps to ns. Another point to note is, when the 
control is turned on or off, it affects not just the fields in directions (1,4) and (2,5). Every 
time the control changes, the field in the third direction pair (3,6) grows and then decays 
to zero. Concurrent to this movement, the intensities in (1,4) and (2,5) undergo a slower 
change, forming a shoulder on each of the time traces. These results are qualitatively 
similar to those in Ref. El but quantitatively dramatically different, in that the intensity 
ratio of the switched signal to the control beam in the present calculation is about two 
orders of magnitude larger than in Ref. IS This improvement is attained by raising the 
radiative decay rate of the cavity from T = 0.4meV in Ref. [8 to T = 1.5meV here, other 
parameters being the same. 

Fig. [7] shows radial (|k|) distributions of intensity in directions 1 and 2 as functions of 
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FIG. 6: (Color online.) (a)-(c) Reflected signals (Ir), summed over radial k modes, vs time 
in the off-axis directions of the multi-|k| model. The directions of the signals in each panel 
are marked by numerical labels, which are defined in Fig. [T£i. Black solid lines represent the 
signals in directions 1, 2, and 3 and the red dashed lines those in directions 4, 5, and 6. 
Imax = 9.1 x 1CT 3 Wcm~ 2 = 3.6 x 10~ 7 /p Ump denotes the peak control intensity, (d) The 
control beam intensity (Iq) vs time. 



time. At about t = 0.2 ns, instability-generated fields appear in a range of radial shells 
in all six directions. At t ~ 1.7 ns, the distribution is narrowed to a single |k| mode, 
which we label by ho, in each direction.^ The polariton frequency of this mode, including 
the effect of the pump-induced HF and PSF shifts, is approximately in resonance with 
the pump frequency. Afterwards, as shown in Fig. |6j the intensity in direction 1 and 
h = h Q rises to a steady state value while that in direction 2 falls to almost zero. The 
control beam is directed at the mode with momentum k = k/j 0j2 . Fig. [7] shows the 
reversible directional switching. As can be seen, there is no spreading of the intensity 
distribution over the radial modes during the switching to direction 2 or the subsequent 
return switching to direction 1. 

The convergence of the off-axis signals in each direction to a single radial momentum 
mode is an important feature of the multi-|k| model simulation results. This implies that 
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a further reduced model, in which only one mode is included in each hexagonal direction, 
would be adequate in analyzing the competitions among patterns in a hexagonal state 
space and the switching process. We call this simpler model the 'single- hexagon model' 
and will use it in the following sections to carry out a parameter- variation study of the 
pattern competitions and analyze their underlying physical mechanisms. In the full 2D 
simulations, the convergence to one radial momentum mode is not as evident, which 
may be due to the fact that the pump beam is a finite-width distribution in k space in 
the 2D simulations instead of a single mode (at k = 0) as in the present model. Results 
similar to those in Fig. |6]have been presented in Ref. El but were not analyzed in detail. 



2. The ring model: results 



In the ring model, we consider modes described by Eq. (13) with k r = 1.2935 x 



10 4 cm 1 , which is the winning radial momentum extracted from the multi-|k| model 
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FIG. 7: (Color online.) Radial distributions of reflected signal intensity in momentum space 
in directions 1 (a) and 2 (b) as functions of time. The intensity distributions in directions 4 
and 5 are similar to those in 1 and 2 respectively. One radial mode always 'wins it all' at 
steady states. The insets indicate the positions of the winning modes in a transverse plane in 
the far field, (c) The temporal profile of the control beam. In this case, the control is directed 
in azimuthal direction 2 with a momentum magnitude matching that of the winning mode in 
direction 1. 
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(Fig. [?]). All parameters are the same as in the full 2D simulations; in particular no 
spatial anisotropy is given to any mode here. Low-level random seed fluctuations in the 
field are put into the off-axis modes. Instability-induced fields initially grow in all modes 
around the ring, but eventually, as in the full 2D simulations, the system stabilizes into 
a hexagonal pattern, the orientation of which is arbitrary and may change from run to 
run. Therefore, the situation that the symmetric ring pattern is unstable, while each 
'broken symmetry' state of a hexagonal pattern is stable, is retained when one reduces 
the full 2D state space to the ring configuration. 

Combining the simulation results in this section, using the full 2D, multi-|k|, and ring 
models, one can see that an important class of transverse instability-driven patterns 
exist in a hexagonal state space. Moreover, a model including only the pump (k = 
0) mode and six modes on a regular hexagon (the single-hexagon model) would be 
sufficient to analyze the competitions and controlled switching among these patterns. 
This analysis will be carried out in the following sections. Interactions among modes 
in different hexagons will be further studied in Section |VIII It will be seen that the 



insight gained from the single-hexagon model is also helpful in understanding these 
inter-hexagon competitions. 

IV. SINGLE-HEXAGON MODEL 

In this and the following three sections, the pattern dynamics in a hexagonal geometry 
are investigated in detail using the single-hexagon model. The results for our standard 
parameter values are reported in this section, and some parametric dependencies of 
the results are studied in Section |V| Detailed analysis of the physical processes driving 



competition and control is contained in Sections |VI| and VII 

The single-hexagon model includes only one sextuplet of points on the vertices of a 
hexagon plus the origin in k space, as illustrated in Fig. |5](c). Within the labeling scheme 
of the last two models we set the number of hexagons N = 1 and drop the index h. The 
momenta of the off-axis modes are written as kj = k r ei,i = 1, ■ ■ ■ ,6, where, as in the 
ring model, k r is assigned the value for the winning mode in the multi-|k| simulations. 
The field symbols are also simplified to % — 0, • • • , 6, with % = labeling the pump 
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fields. Figs. [8] and [9] show the simulation results using the same material and run-time 
parameters as in the multi-|k| simulation presented in the previous section. As in that 
simulation with the multi-|k| model, modes 1 and 4 are given a slight advantage, and 
a control beam is periodically applied to mode 2. Exciton densities in three off- axis 
= 1, • • • ,3, are plotted as functions of time in Fig. pBl Their behaviors 



directions, |pi| 2 ,i = 1, • 
are clearly similar to those of the reflected light fields obtained within the multi-|k| 
model (Fig. [6]). In the winning off-axis direction, the steady state exciton density is of 
the order of 10 8 cm -2 , or about 0.01 of the pump induced exciton density. While the 
rise time of pi in each cycle is about the same as in the multi-|k| model results, the time 
for complete switching, when the control is turned on, from direction 1 to 2, is roughly 
200 ps longer here. As explained in the next subsection, a control intensity threshold for 
complete switching exists, and the longer switching time in the single-hexagon model 
indicates that this threshold is higher here than in the multi-|k| model. 



Since the polariton scatterings represented in Eqs. (14) and (15) are coherent pro- 
cesses, the density transfer among modes depends critically on the relative phases of the 
fields involved. We display the phases of off-axis exciton fields relative to the pump's 
phase in Fig. [9j Their behaviors and their roles in the dynamics of pattern selection will 
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FIG. 8: (Color online.) Exciton density for three off-axis directions in the single-hexagon 
model: \pi\ 2 - black solid line, \p2\ 2 - red dashed line, \ps\ 2 - blue dotted line. Note that the 
dotted line, |ps| 2 , overlaps with the other lines over much of its range. The control beam, with 
peak intensity 3.6 x 10 _7 I pump is first switched on at t = 2 ns in this run. 
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be more thoroughly discussed in Sections VI and VII below. The importance of phases 
was also investigated in other instability-driven structures, e.g. in Ref. H6] 

The exciton density in the pump direction (|po| 2 , which is not shown) reaches its 
peak value of about 10 10 cm~ 2 at around t = 30 ps and stays constant except for a 
slight variation when the control is on. The phase of the pump exciton field also stays 
essentially constant at roughly 0.157T after t = 30 ps. 

The pattern-formation time starting from the initial pump-only state varies with 
the models used. In general, this initial formation time is longer when there are more 
competing modes in the model, being shortest in the single-hexagon model and longest in 
the full 2D simulations. The switching time when the control beam is applied, however, 
is roughly the same in all models. 
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FIG. 9: (Color online.) Phases of the exciton fields in the off-axis directions for the reference 
case within the single-hexagon model, (a) Phases in modes 1-3, and (b) phases in modes 4-6, 
relative to the pump-induced exciton phase, Scfrj = <pj — 4>o,j = 1 — 6. The modes are labeled 
in Fig. [5]3. Note that the symmetry between <j>j and 4>j+3 is broken after the introduction of 
control beam at about t = 2 ns. 



29 



V. PATTERN AND TIME SCALE VARIATIONS WITH PARAMETERS 



In the previous two sections, we have discussed a particular scenario (with asymmetry 
applied) where a 'two-spot' pattern is selected by the system's dynamics and then the 
pattern is rotated by a control beam. In this section, we investigate the range of stable 
patterns supported by the hexagonal state space and how the selection of patterns as well 
as the transition time between patterns depend on the system and control parameters. 
We will refer to the set of results presented in the previous section as the 'reference case'. 

A. Varying the anisotropy in the cavity mode energy 

In the reference case, we introduced an anisotropy in the cavity mode energy cj£ by 
down-shifting it by 0.12 meV in modes (1,4). We examine here the effects of changing 
the size of this anisotropy, which we denote by 5. The reference value of 0.12 meV 

4 

3 

o 2 

00 

O 

T 1 



0.0 0.4 0.8 1.2 

5/5 

FIG. 10: (Color online.) Steady state exciton densities in directions 1 and 2 vs the amount of 
anisotropy - downward shift of the cavity mode energy in direction 1 (5). 5q = 0.12 meV is the 
shift in the reference case (see text). |^>3 1 2 (not shown) is equal to \p2\ 2 - For < 5 < 5p ~ 0.285o, 
only the hexagonal pattern is stable. For 5 > 5 U ~ 0.66<5o, only the two-spot (in modes 1 and 
4) pattern is stable. Both patterns are stable in the range between 5g and 6 U , where their 
adiabatic evolutions show hysteresis behaviors: the solid (dashed) lines trace the evolution of 
the densities as 5 increases (decreases) from small (large) values. 
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is denoted by Sq. Figure 10 shows the exciton densities in directions 1 and 2 in the 
steady state (s) as a function of 6/ So. The density in direction 3 is equal to that in 
direction 2 in each case. For S = 0, where there is complete symmetry among the six 
directions in the equations, we have found only the symmetric steady state solution, 
with \pi\ = \p 2 \ = \ps\, to be stable. 52 The hexagon remains the only stable pattern for 
values of S between and a lower critical value Se ~ 0.285o, with the steady state value of 
\pi \ steadily growing while that of \p2\ = \ps\ is falling. Between Si and an upper critical 
value S u ~ 0. 665q, both the hexagon and the two-spot pattern, in directions (1,4), are 



stable. Above S u only the two-spot pattern is stable. In Fig. [TOj the solid lines trace 
the evolution of the steady state densities if S is slowly increased from zero: they vary 
smoothly until S crosses S u , when \px\ jumps up and |^> 2 1 (= \Ps\) drops abruptly to zero. 
Starting from large S, when we decrease S, the system follows the dashed lines inside the 
bistable range. The total off-axis density, summed over all six directions, increases with 
S over the displayed range (not shown). For example, the total exciton density in the 
symmetric case (S = 0) amounts to about 50% of that in the reference case (5 = Sq). 
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FIG. 11: (Color online.) Exciton density in direction 1 (preferred direction) vs time for various 
amounts of downward shift in the cavity mode energy (5) calculated in the single- hexagon 
model. The value of S/5q for each curve is given in the inset. For 5 below the critical value 



for hexagon-to-two-spot transition (see Fig. 10), attains steady state within 1 ns. For 5 
immediately above this critical value, e.g. the curve for 8 /So = 0.67, \p\\ passes through a long 
intermediate step before rising to its steady state level. 
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In Fig. 11, we compare the time evolution of \pi\ 2 for various values of S during 
the transition from the initial vacuum state to the stable steady state. For this initial 
state, the system goes asymptotically to the hexagon steady state for values of 5 in the 
bistable range. In each case, the pump intensity becomes practically constant at around 
t = 30 ps. The off-axis density becomes visible on the scale of the final steady state level 
at about 0.6 ns. After this point, for 5 < 5 U , \pi\ 2 takes a relatively short time - several 
tens of ps - to reach its steady state value. Beyond the hexagon-two-spot transition, 
however, the time trace of |pi| 2 acquires an intermediate step, during which its change 
slows down, before its final rise to a constant value. The length of this intermediate 
step appears to increase without limit as 5 approaches S u from above. Similar to the 
reference case, shown in Figs. ^ and for each 5 > 5 U , |p2^)| 2 and |p3(i)| 2 stay non- 
zero during the time when |pi(t)| 2 is passing through its intermediate step. Comparing 
the curves at 5 = 0.65<5o and 5 = 0.67<5o suggests that the intermediate step is a 'residual' 
effect of the hexagonal pattern. Though no longer a stable steady state, it still slows 
down the system's dynamics when the system passes through its vicinity. In cases not 
close to the hexagon-two-spot transition, one lower bound to the formation time of the 
stable pattern is the ps time scale of the electron-hole Hamiltonian. In most cases, 
the actual time scale is one to two orders of magnitude longer that this lower bound 
depending on the level of fluctuations and the balance between gain and (both radiative 
and non-radiative) loss rates. 

When the asymmetry 5 is set below 0, directions 1 and 4 are given a disadvantage, 
and the stable pattern is a hexagon with \p 2 \ 2 = |jo 3 1 2 > \pi\ 2 . A transition to a four-spot 



pattern would not occur since the latter is not a steady state solution to Eq. (15) (see 



Section VI below). 

The switching behavior of the hexagonal pattern for 5 < 5 U is similar to that in the 
reference case when the same control beam is applied. In the symmetric case (5 = 0, 
for example), when the control beam is turned on in direction 2, the hexagon switches 
to the two-spot pattern with the density in directions 2 and 5 being the same as in 
the reference case. When the control beam is then turned off, the system goes back to 
the symmetric solution. The times needed for switching between the two patterns are 
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slightly shorter than those in the reference case. 



B. Varying the control beam intensity 



In the reference control beam with peak intensity = 9.1 x 10 3 Wcm 2 = 

3.6 x lCF 7 / pump is applied in direction 2. We examine here the effects of varying this 
control intensity, retaining for all other parameters their values as in the reference case. 
We denote the control intensity by J max and its value in the reference case by i^ax- 



Fig. 12 shows the steady state exciton densities when the system is under irradiation 
by the control beam as functions of / max - It can be seen that the pattern undergoes two 
abrupt changes. When the control is off (I m3x = 0), the stable pattern has two spots 



in directions (1,4). As J max increases from zero, |pi| 2 (=|p4| 2 ) falls while \p2\ 



and | jo 3 1 2 (=|p6| 2 ) rise gradually. As J max reaches a critical value / CT (i), which is between 
and /max' the hexagon with two bright spots and four dim spots changes abruptly to 
a hexagon with comparable brightness in all spots. Thereafter, \pi\ 2 continues to fall 
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FIG. 12: (Color online.) Steady state exciton densities vs the peak intensity of a control beam 
introduced in direction 2. The control intensity (7 max ) is expressed in units of that in the 
reference case, denoted here by /max- As W increases, the pattern goes from two spots (in 
directions 1 and 4) through a hexagon to two rotated spots (in directions 2 and 5). A bistable 
region and attendant hysteresis behavior (not shown in the figure) are present at each abrupt 
transition. 
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while |p 2 | 2 continues to grow until J max reaches a second critical value J C r(2), which is 
0.86/^ ax , where |pi| 2 and | j>3 1 2 drop to zero and the pattern changes abruptly to two 
spots in directions (2,5). Hence I C r(2) is the lowest control intensity that can completely 
switch the two spots from directions (1,4) to (2,5). Though not shown in Fig. 12, a 



bistable region and hysteresis behavior, similar to those in Fig. 10, are present at each 
abrupt transition. 

The time traces of the densities in the reference case, in Figs. [6]and|8j show interme- 
diate steps also during the switching process. As in the initial pattern formation stage, 
the intermediate step here can be interpreted as a 'transient trapping' of the system in 



the vicinity of the hexagonal state that loses its stability at I D 



'cr(2) 



. The duration 



of this intermediate step also becomes longer as we set the control intensity closer to 
the critical value I C r(2)- For J max sufficiently close to I C r(2), the time the system spends 
in this transient trapping dominates the time taken by the two-spot pattern to switch 



directions. Illustrating this, we plot the switching time as a function of I max in Fig. [13 
Again, there does not appear to be an upper bound to this duration. Away from 7 cr (2) 
when the system is no longer affected by the hexagonal state, the switching time is 
around 1 ns. 

If the asymmetry 5 is decreased, the critical control intensity for rotating the two-spot 
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FIG. 13: Time taken to switch between a two-spot pattern in directions (1, 4) to one in 
directions (2, 5) vs the peak intensity of a control beam introduced in direction 2. The control 
intensity (7 max ) is expressed in units of that in the reference case, denoted here by Im ax . 
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pattern I C r(2) is also lowered. However as there is a minimum value for 5, namely S CI for 
the two-point solution to be stable, a lower bound to the switching intensity exists for a 
fixed set of material parameters and pump intensity. 



VI. DYNAMICAL ANALYSIS 



In the following two sections, we seek to gain insight into the physical mechanisms 
underlying the phenomena of pattern formation and switching displayed in the previous 
sections. For this purpose we analyze the interplay between the various wave mixing 
(or coherent scattering) processes contributing to the polariton dynamics. To keep the 
discussion most transparent, we will work within the single-hexagon model. 

It is convenient to break the complex fields in the rotating frame down to their 
magnitudes and phases and consider the equations of motion for these components: 

Pj = p je -^ Ej = Eje-ite-^ (16) 

for j = 0, • • •, 6. Here £j is defined to be the phase difference between the exciton field 



Pj and the cavity photon field Ej. Substituting Eq. (16) into Eqs. (14) and (15), we get 
for the off-axis fields (j = 1, • • •, 6), 



dt 
dcpj 
~dt 



— lm[ihpjp*] = L Pj + Q p . + C p . 
npj 

j^Re[ihpjP*) = j_ [L+. + Q+. + CVJ 



(17) 
(18) 
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where 



and 



hQ 



pj 



nc, 



HL p . = V RF p j+3 Po sin ($°'° +3 ) + 2Anp E p j+ 3 sin ($JJ +3)So 



o.o 



SlEj sin (£_,-)(! - 2Apl) + 2AVtp Q E p j sin (f 



IxPj 



2p 



':-■:> V m popj-i sin + ifi-BoPy-i sin 

yU>/ A) /; y ; sin !>I>^ 



+ 2 Pj+2 [VhfMj+i sin ( $ i,i+2 



+ AQE p J+1 sin ($J^ ) + AQp E J+1 sin (Sjft^; 



+2p 



vkFPi+ift-i sm ci-;;/-' ') + in^+^-i sin ($^- ] 



Sj+1 



+ Alpi+i^-i sin (^J^ 1 ) 
2VbPj+3 sin ( 'I' ;.}';; 2 ) + Pj+2pj-i sin (fcJ-J^ -1 

+2lOp, +3 £ p^ +3 sin($^ i6+3 ) 

sin (£,) ^ + 2AVtpj ^ Pi E i sin (&) 



ftL^ = [ejjj - ftcdp + Vb-Po] Pi + ^HFPi+3Po cos 
+2AnpoE p j+3 cos (^Jj +3j& ) 



0,0 \ 



VLEj cos (^) 1 - 2Ap z Q + 2v4fip ^o cos (£ )Pj 



where ^J'" 1 



'i + 01 



ft<3 Pj - [sin i — y cos] 

HC Pj [sin h-> cos] + VuFpj (2 ^ p 2 + 



(19) 



(20) 



(21) 



and $ 



l,m 



(22) 
(23) 
(24) 

In these equations, 



each sum over i ranges from 1 to 6. In Eqs. (23) and (24), a function followed by the 



symbol [sin y cos] means the same function with all sines replaced by cosines in its 
expression. We have grouped the contributing terms according to their order in the off- 
axis fields. As the symbols imply, L p ., Q Pj , and C p . contain respectively first (linear), 
second (quadratic), and third-order (cubic) terms in pj or Ej, j — 1, • • -,6. The same 
applies to , Q^. , and C ( p j . The state with polariton density only in the pump mode 
is treated as the zeroth-order state. The pump fields po and Eq experience very small 
(fractional) changes in time and will be taken approximately as constant in this analysis. 
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As mentioned before, all exciton interactions and PSF terms can be visualized as 
either a four-wave mixing process with wave-vector matching or a polariton-polariton 



scattering process with momentum conservation. For each scattering process in Eqs. ( 14 ) 



and (15), the exciton field on the left- hand- side and the conjugated field on the right- 
hand-side are in the outgoing modes while the other two fields on the right-hand-side are 
in the incoming modes. Of course these scatterings are coherent, not of the Boltzmann, 
mass-action type: the direction of density transfer depends on the interference between 



the scattered wave and the existing field in each outgoing channel. In Eqs. (17)-(24), 
the phases of the incoming fields are subtracted from those of the outgoing fields inside 
the sine or cosine functions. Since the density equations contain the sine functions, the 
term is a gain term for the mode in question if the combined phases of the outgoing 
fields lead (by less than it) those of the incoming fields. Otherwise, it is a loss term. 
For the phase's rate of change, with the cosine functions, a process tends to increase the 
phase in question if the phase difference between the outgoing fields and the incoming 
fields is less than ir/2. 



The first two terms contributing to L p . in Eq. (19) are the primary instability-driving 
(gain) processes, in which two pump polaritons scatter off each other into two opposite 
off-axis directions. An example is displayed pictorially in Fig. |l4|a). The next term is 
the exciton-photon coupling in an off-axis mode, partially blocked by the pump exciton 
density. In the calculations, the angle £ 3 - is always positive, corresponding to an expected 
net energy flow from the exciton field pj to the photon field Ej. Hence, to the exciton 
field, this term represents the radiative loss from mode k,. The next term describes 
the forward scattering between a pump photon and an off-axis exciton, or equivalently, 
a self-wave-mixing process in which a pump photon scatters off a density grating, set 
up by the pump and the off-axis exciton fields, into the off-axis direction. In the pump 
mode, energy flows on balance from the cavity photon field to the exciton field, resulting 
in a negative phase-lag angle £ - Hence, this scattering actually results in a density 
transfer from the off-axis k, mode to the pump mode. The last term — IxPj represents 
all non-radiative losses. The balance between gain and loss terms determines the sign 
and magnitude of L Pj . 
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In a quadratic process, i.e. one that contributes to Q Pj and Q^, one off-axis polariton 
(from mode kj) and a pump polariton scatter into the modes kj_i and k J+ i and vice 



versa. These processes are also illustrated in Fig. 14 the first three terms in Eq. (20) 



are represented by Fig. 14 b), the next three terms by Fig. 14 c), and the last three 



terms by Fig. |14|(d). We note that the momentum conservation condition that gives rise 
to the quadratic terms is satisfied only for modes arranged geometrically in a hexagon. 
Among the cubic terms, the first three contributing to C p . and . represent two counter- 
directed polaritons on the hexagon scatter into a different pair of opposite directions. An 
example is drawn in Fig. [l4|(e) . The next two terms are, respectively, the Pauli blocking 
reduction of the (linear) exciton-photon coupling by the exciton density in all off-axis 
modes and a self- wave-mixing process in which a photon in mode i (i ^ j) scatters off 



(a) 






FIG. 14: (Color online.) Pictorial representations of polariton scattering processes contributing 
to the equation of motion of p%. The processes are (a) linear, (b), (c), (d) quadratic, and (e) 
cubic in the off-axis fields. In each diagram, the incoming modes and outgoing modes in 
the scattering are labeled by red, solid circles and blue, open circles respectively. The arrows 
represent the modes' momenta, dashed (solid) for the incoming (outgoing) modes. For example 
in (b), two polaritons from directions (pump) and 6 scatter into directions 1 and 5. 
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density grating, set up by pi and pj, into mode j. Since < £j < 7r/2 for i = 1, • • •, 6, 



these two terms are positive in both Eqs. (21) and (24). The last term contributing to 



Cfj,. in Eq. (24) is a (blue) shift in frequency for pj. In contrast to the quadratic terms 
discussed above, these cubic processes do not require a hexagon mode configuration to 
be operative. 



We finally note a phase degeneracy in the solution to Eqs. (17)-(24). Suppose 
(pj(t), 4>j(t)), j = 1, • • • ,6 is a solution to these equations. Then one can verify that 
(pj(t), 4>j(t) + 5<pj) is also a solution provided the six (50/s are time-independent and 
satisfy the following constraints: 

5(p j + 5(j)j + 3 = 271171 (25) 
S(j)j + i + <50j_i — 5<pj = 27im (26) 

where m is an integer, and the subscript is again counted cyclically through 1 to 6. 
These two constraints restrict the number of undetermined phases to two. This applies 
in particular to steady state solutions. Thus there exist an infinite number of steady 



state solutions to Eqs. (17)-(24) which can be grouped into a finite number of classes. 
The solutions in each class have the same set of magnitudes and are parameterized by 
two free phases. We see an example of this phase freedom in Fig. [9j when we compare 
the phases of, say, p\ and p^ during the initial steady state with their phases when the 
system returns to the (1,4) two-dot pattern after the first on/off switching cycle. The 
individual steady state phases are different in the two time periods, while their sum stays 
the same. Therefore, strictly speaking, the system 'returns' after the on/off switching 
cycle to a different steady state that has the same set of magnitudes as before. Put 
another way, in a hypothetical experiment, the measured relative phase between E\ and 
Eq may be affected by uncontrollable factors, but the measured average phase of Ei and 
E4 relative to E'o's phase can be meaningfully compared to our theory's prediction. 



VII. INTERPLAY OF WAVE MIXING PROCESSES 



We now discuss the contributions and interplay of the various processes in effect- 
ing pattern selection and switching in the reference case (Section IV). It is clear from 
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Eqs. ( 19 )-(21 ) that to understand how the exciton densities evolve, one needs infor- 
mation on the relative phases of the polariton fields involved in the scatterings. The 
calculated magnitudes and phases of the off-axis pj for the reference case were plotted 



in Figs. [8] and [9j Figure 15 shows the contributions of the three groups of terms - linear, 
quadratic, and cubic - to the rates of change of the exciton densities in directions 1 



and 2, and Fig. 16 shows these contributions to the evolution of the respective exciton 
phases. Initially, the linear terms dominate, and in each off-axis mode, the exciton phase 
4>j is quickly locked into a value given by the solution of d(j)j/dt ~ L^/pj = 0, where 



is given by Eq. (22). Actually, as explained above, this condition fixes only the value of 
the sum <pj + J+3 but not either phase angle individually. In the calculations, however, 
all the exciton phases are initially set at 0, and since L^/pj « 3 /pj +3 , we also have 
4>j ~ 4>j+3 in the linear regime. The radial |k| value of the off-axis modes having been 



chosen optimal, we have (cf Eq. (19)) 



d>J° +3 = (j>j + i+3 - 20o ~ tt/2. (27) 

In Fig. [9j this phase lead is about 0.67T for each of the three off-axis mode pairs. Once 
their phases are locked, the off-axis exciton densities grow exponentially from fluctuation 
levels. We are considering pump intensities close to the phase-conjugate instability 
threshold. This condition implies a substantial cancelation between the gain and loss 
terms in L p .. In the reference case calculations, the exponential off-axis growth rate in 
the first 0.6 ns is about 0.01 ps -1 , which is much slower than the typical inverse time 
scales given by the exciton Hamiltonian (e.g. ^x/h = 0.67 ps _1 ). During this exponential 
growth phase, the p/s are actually too small to be visible on the scale in Fig. [8j It is 
interesting to note that, although the asymmetry-favored p\ (and p±) eventually wins out 
at steady state, p 2 actually grows slightly faster than p 1 during the initial 0.6 ns. Thus 
the asymmetry confers its advantage through the quadratic and/or cubic processes. Put 
another way, the selection of (pi, p^) as the winning directions by the asymmetry works 
through the interactions among excitons in all six off-axis directions. 

The linear growth regime lasts approximately 0.5 ns. As the p/s grow, the quadratic 
and cubic contributions to dpj/dt come into play. The onset of the quadratic terms has 



two effects. First, their direct contribution to the magnitude equation, Q p in Eqs. (17) 
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and (20), is positive, as shown in Fig. 15 This can be understood as follows. As 



explained in Section |VI[ whether a scattering term is a gain or loss process depends on 
the phase sum of the outgoing exciton-polariton fields relative to that of the incoming 
fields. In each term contributing to Q p , one of the four exciton fields is in the pump 



mode while the other three are off-axis. As we have seen above in Eq. (27), the linear 



instability sets up the off-axis fields with phases leading the pump's by a large margin. 



In the first six terms of Eq. (20), represented by Fig. |l4[b) and (c), po is among the 
incoming fields, and thus the relative phase angles are typically positive (between and 
7r). These six terms are therefore gain terms. The photon/exciton phase shift in each 
mode, is very small and its effect can be ignored in this argument. In contrast, the 



last three terms in Eq. (20 ), represented by Fig. 14 cl), yield losses because p is outgoing 
in them. Since the three triplets of terms are roughly equal in magnitude, and the gains 
outnumber the losses, the net result for Q p . is typically positive. 

The second effect of the quadratic terms acts through their contributions to d<pj/dt 



in Eq. (18). Again from the linear phase lead, Eq. (27), one can deduce that the relative 



phases in the terms contributing to Qs in Eq. (23) typically lie, during the onset of the 



quadratic terms, between — 7r/2 and tt/2, Q^. is therefore positive, as shown in Fig. 
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The resulting increase in </>,- exerts a negative feedback effect on the linear scattering 
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FIG. 15: (Color online.) Contributions of the linear, quadratic, and cubic terms to the rate of 
change of exciton densities in directions 1 (a) and 2 (b) in the reference case. 
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°'° +3 ) is a little 



process in L p : since the relative phase in the linear gain terms sin($ 
over 7r/2 during the initial growth period, pushing cf)j (and 0^+3) up reduces the gain. 
And since the loss terms, e.g. •y x pj, are not affected, L Pj suffers a net decrease and 



becomes negative, as can be seen in Fig. 15. In fact, the upsurge of pi/4 to their winning 
(two-spot) steady state is associated with a slower rise in 0i/ 4 than in the other modes' 



phases. Initially locked at zero, La. is also decreased by the growing 3 negative 



JJ+3 



values, as shown in Fig. 16, thus offsetting the effects of Qa. to some degree. 



The contributions of the third-order scattering processes to the exciton density, i.e. 



the first three terms in Eq. (21 ), depend on the relative phases between the off-axis fields. 
When pi is winning, since 0i < 02 = 03, inspecting Eq. (21) shows that these terms are 
negative for C Pl and positive for C P2 and C P3 , or they tend to equalize the densities in 
the six modes. The remaining two terms in C p . are Pauli blocking terms that reduce the 
flow of density from pj to Ej and so are gain terms to mode j. In the two-spot steady 
state with modes 1 and 4, the only surviving contribution to C Pl is 2A£lEip\ sin (£1). 



Similar to Qa., all terms in Ca. in Eq. (24) are positive, thus again helping to lower the 



linear amplitude term L Pl . The first set of terms in Eq. (24) again represent coherent 
scatterings and Pauli blocking. The last term is a blue shift of the exciton energy in 
mode j due to interactions with off-axis excitons. 
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FIG. 16: (Color online.) Contributions of the linear, quadratic, and cubic terms to the rate of 
change of exciton phases in directions 1 (a) and 2 (b) in the reference case. 
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The above analysis underlines the fact that pattern competition in even the single- 
hexagon model is the result of rather complex, activating or inhibiting, feedbacks among 
the various physical processes. Not being simply additive, the effects of the processes are 
not easy to disentangle cleanly. Nevertheless, it is instructive to consider an artificial 
reduction of the single-hexagon model in which the quadratic processes are disabled. 
In this case, the most basic function of the cubic processes is to saturate the polariton 
density in the off-axis modes, which the cubic processes achieve mostly through raising 
the phases of p/s, as explained above. Setting the quadratic terms (last three lines in 
the single-hexagon of Eq. (15), or equivalently, Q Pj and in Eqs. (17) and (18)) to 
zero, we have repeated the simulations with the same parameters as in the runs shown 



in Section [Vj Fig. 17 shows the exciton densities in modes 1 to 3 for an isotropic setting, 
i.e. no advantage is given to modes 1 and 4 (5 = 0). Unlike the 'physical' case, with the 
quadratic terms operative, where the symmetric hexagon pattern is stable, the stable 
steady state in Fig.[l7]is a two-spot pattern - the six-fold symmetry of the single- hexagon 
model is spontaneously broken. Which pair of modes win at the expense of the others is 
not determined by the model's physics and changes from run to run. When a sufficiently 
strong control beam is applied to mode 2, the pattern switches over to modes 2 and 5, 
but stays in this pair when the control is turned off. This is consistent with the fact 
that, in the isotropic setting, the two-spot pattern at modes 2 and 5 is the stable state 
closest to the steady state maintained by the control beam. 

It is useful here to draw some insight from previous analyses of other multiple-mode 
competition modelsS^SSl that are based on linear growth and cubic saturation. Consider 
Lamb's two-mode laser modelj^where each mode exerts a cubic saturation effect on both 
itself and the other mode. For this model, a symmetric steady state, with equal intensity 
in the two modes, is stable if self- saturation is more efficient than cross-saturation, 
while two asymmetric states, each with intensity in only one mode, are stable if cross- 
saturation is more efficient. By analogy, we also subdivide the contributions to C^. 
and C Pj into those due to densities in modes in j and j + 3 and those due to densities 
in other modes, and refer to the two classes as 'self-saturating' and 'cross-saturating' 
contributions respectively. The analogy is not exact because the effects of the two classes 
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on dpj/dt are not simply additive, as their counterparts in Lamb's modeP^are, but these 
arguments are still applicable qualitatively, and our numerical results indicate that the 
cross-saturation cubic processes in our model polariton system are more efficient than 
the self-saturation processes. We have already seen that turning the quadratic processes 
back on in the isotropic setting destabilizes the two-spot pattern while restoring stability 
to the symmetric hexagonal pattern. This indicates that the enhancement brought by 
Q p , to dpj/dt outweighs the suppression caused by the positive Q^. in raising the phase 
of pj. Thus the net effect of the quadratic processes can be classified as 'cross-activating'. 
The introduction of anisotropy (5 > 0) modifies the balance of the above activating and 
saturating effects and, when strong enough, destabilizes the hexagonal pattern. 

In the above, we have analyzed the physical processes driving transverse pattern 
selection and control in a quantum well microcavity and attempted to pursue a semi- 
phenomenological analogy between the pattern dynamics in this system and other mode- 
(or pattern-) competition systems. We will explore this analogy further in a future pub- 
lication. 




Time [ns] 

FIG. 17: (Color online.) Exciton densities in directions 1 to 3 for the reference case but with 
the quadratic terms disabled and no asymmetry in the cavity dispersion relation. The first 
on/off cycle of the control beam starts at t = 7ns. 
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VIII. COMPETITIONS AMONG HEXAGONS 



With the single-hexagon model, we have analyzed the switching between patterns 
that are subsets of the same (regular) hexagon. In this section, we investigate switching 
between patterns that reside in different hexagons, as displayed in simulations using 



the multi-|k| and ring models. One can see by inspecting Eqs. (15) that the quadratic 



processes (last three lines of Eqs. (15)) take place only among modes situated on the 



same hexagon. Modes on different hexagons compete with each other only through the 
cubic processes. From the intuition that we formed in the previous section, that cubic 
cross-saturation in our system is more efficient than self-saturation, we expect patterns 
that break the symmetry among hexagons to be favored. We have already seen that 
this is consistent with the simulation results in the 2D r-space model and the ring model 
in Section |III[ In both simulations, the setup is isotropic, but a hexagon pattern with 
arbitrary orientation is stable while the ring pattern is not. Within the hexagon, the 
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FIG. 18: (Color online.) (a) Time evolution of the angular distribution of the reflected signal 
intensity in momentum space, calculated in the ring model. Although the setup is perfectly 
isotropic, the ring pattern is unstable while the hexagons are stable. The orientation of the 
hexagonal pattern obtained in each run is arbitrary. A control beam is introduced at t = 7ns, 
directed at a mode outside of the initial hexagon, switching the pattern to a two-spot one. 
When the control beam is switched off at t = 9.5ns, the signals form a hexagon in the orientation 
defined by the control beam, (b) Temporal profile of the control beam intensity. 
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cross-activating quadratic processes stabilize the symmetric hexagonal pattern at the 



expense of the two-spot solutions. In Fig. [18] we show some extended results of the ring 
model simulation. The figure shows the time evolution of the angular distribution of 
intensity. Starting with only the pump mode, the system forms a hexagon at about 
t ~ 6.5ns, and the exciton density in the winning hexagon is 5.6 x 10 7 cm -2 , which 
is the same as the steady state density in the symmetric multi-|k| model (not shown). 
A control beam with the same intensity as the one in the reference case is applied to 
another mode on the ring. Just as in the reference case, the control stabilizes a two-spot 
pattern with one spot in the mode it shines on. When the control beam is subsequently 
turned off, the off-axis pattern reverts to a hexagon which is oriented to contain the 
previous two-spot pattern as a subset. This is again consistent with the fact that each 
individual hexagon pattern is stable, and when the two-spot pattern is rendered unstable 
by the control beam's removal, the system settles into the nearest stable pattern. 

A comparable situation is present in the multi-|k| simulations. We performed a 
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FIG. 19: (Color online.) Time evolution of the radial distributions of reflected signal intensity 
in momentum space in directions 1 (a) and 2 (b), calculated in the multi-|k| model. One radial 
mode always 'wins it all' at steady state. The inset in each panel indicates the positions of the 
winning modes in a transverse plane in the far field, (c) Temporal profile of the control beam. 
In this case, the control is directed in azimuthal direction 2 with a momentum magnitude 
smaller than that of the winning mode in direction 1. 
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simulation similar to the one shown in Figs. [6] and [7] except that the control beam 
is applied to a mode on a different hexagon (i.e., having a different radial momentum 



|k|) from the initial two-spot pattern. Recall the notations in Section III with the 
anisotropy-induced preference, the initial two-spot pattern is in modes k ft0jl and k/i 0j4 , 
and the control beam is applied to mode k^ 0)2 . Here the control is directed at mode k^ )2 
with h' being a close neighbor of h . Fig. 19 shows the time evolution of the exciton 
density in modes over a range of h, or |k|, values in directions 1 and 2. One can see that 
when the control is on, the pattern follows it to the mode pair k^ 2 and k^ 5 . When 
the control is turned off, the pattern switches back to directions 1 and 4 but stays in 
the radial momentum mode set by the control, i.e. to the modes k/^,2 = 1,4 instead 
of the original k hQti ,i = 1,4. In subsequent cycles of switching, the pattern stays in 
the hexagon defined by h' . These behaviors again indicate that, in the absence of the 
control, there exist a group of stable steady states, in each one of which the signal is 
concentrated in a mode with k = \th,i,i = 1,4 and h is within a certain range about 
ho. Which state the system evolves into depends on the system's history, e.g. which 
mode the signal is in when the control is turned off. In contrast to the situation in the 
ring model, however, there is no symmetry under which these modes are equivalent: for 
example, they have different single-polariton energies. 



IX. SUMMARY AND OUTLOOK 



In this paper, we have studied the nonlinear polariton dynamics of a laser-pumped 
quantum well microcavity as a spontaneous pattern-forming system. Using a microscopic 
theory of exciton-exciton interactions and GaAs parameters, we have examined the 
formation, selection, and control of the various transverse polariton density modulation 
patterns. We have analyzed in detail the polariton scattering processes that govern 
the competitions among patterns, obtaining useful insights into the crucial role played 
by the (relative) phases of the polariton amplitudes in driving the gains and losses 
of the densities in various k modes. We have also performed a more comprehensive 
investigation of a previously proposed low-intensity, all-optical switching scheme^ that 
exploits these patterns. 
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One outcome of our analysis is a qualitative characterization of the effects of the var- 
ious scattering processes in terms of notions used in discussing pattern formation and 
competition in other, usually classical, physical systems. These notions are based only 
on the polariton densities, despite the fact that the competitions among polariton am- 
plitudes acts to an important extent through their phases. This has helped us form an 
intuitive picture of the pattern dynamics in our polariton system. We have pursued this 
analogy with other pattern-forming systems further and have constructed a three-state 
mode-competition model which can be related to our models here and whose dynamics 
mimic almost all the qualitative behaviors displayed so far by the solutions of the single- 
hexagon model. Because of its simplicity, this mode-competition model is amenable to 
semi-analytic treatment, and we have obtained a complete classification of the steady 
state solutions in the relevant region of parameter space. This 'global' perspective in- 
troduces a useful organizational framework to discuss the qualitative dynamics of our 
system. This work will be reported in a future publication. 
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Appendix A: Coupling between the cavity modes and macroscopic fields 

In this appendix, we give more details on our model for the coupling between the 
macroscopic light field and the cavity modes. The relations Eqs. ([5]) - (J6j) between input 
and output fields will be derived, and the choice of Eq. ^ will be given support by 
showing it leads to the interpretation of the cavity field magnitude squared \E^\ 2 as 
a photon density. In our model, the cavity modes are taken to be the spatial Fourier 
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modes of an oscillator field separate from the macroscopic light field and confined to an 
infinitely thin plane. The QW is assumed to be embedded inside this planar cavity, and 
the oscillator field and the exciton field are coupled to form the polariton field. The 
(linear) coupling between the cavity oscillator field and the outside radiation field is 
modeled like the coupling between the latter and a material oscillator field. Since all 
macroscopic light modes involved are near-normal (deviations of their wave vectors from 
the normal direction are typically less than 10 degrees), for simplicity, the complications 
due to the finite polar angles of the off-axis modes are ignored. More realistic modeling 
of the microcavity and its coupling to macroscopic light field would lead to quantitative 
corrections to our results but are not expected to change any conclusions in this paper. 

The propagation of the macroscopic light field is governed by the Maxwell equations 
(with no external 'free' charges or currents). The medium outside the cavity is taken 
as a dielectric with a real refractive index n s = and magnetic permeability fi s = 1. 
The oscillator field acts as a local contribution to the polarization density, yielding the 
following constitutive equations: 

V(r, t) = e n 2 s £(r, t) - M c E cav (x, y, t)S(z) (Al) 
H(r,t) = —B(r,t) (A2) 

Here the functions £, V, B, and H, denote the positive-frequency parts of the respective 
fields so that, e.g., the (real-valued) electric field is given by 2Re£. E cav (x,y,t) is the 
amplitude of the resonance mode inside the cavity whose equation of motion is given 
below (and in Ref. [8]), and t c is a parameter characterizing the coupling between the 



light field and the oscillator mode. Substituting Eqs. (Al) and (A2) into the Maxwell 
equations and eliminating B in the usual way, we obtain a wave equation for the radiation 
electric field £: 

— V (V • £ ) + n 2 s e ^£ - -V 2 £ = M c 5{z)^E c&y (A3) 

For near-normal fields, V • £ ~ since V • £ is exactly zero at any point outside the 
cavity, and, across the cavity, it is proportional to the difference in E z which is small. 
We take the approximation V ■ £ = in the above equation, and replacing eo/^o by 1/c 2 
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we obtain 

n 2 d 2 fit d 2 

We expand the fields in their spatial Fourier modes in the cavity's plane, i.e., in the 
(x, y) directions: 



^M) = ^e^ + ^4(2,f) (A5) 
L k 

E cav (x,y,t) = ^^e^ +fc ^E cav , k (t) (A6) 



where k = (k x , k y ) and C 2 is the planar area of the cavity. Substituting into Eq. (A4) 
we have for each Fourier component 

n 2 d 2 ( d 2 \ fit d 2 

where k 2 = k 2 + k 2 . The (positive frequency part of the) general solution of this equation 

at any point outside of the cavity, z ^ 0, is given by: 

f°° di ) r ~ i 

£ k (z,t) = — e i( ^ (fc ' w)2 - wt) 5 k +(w) -e 4( - fcz(fc ^- wt) 4-(w) (A8) 

Vn 27T L ' - 



where k z (k,u) = +\j'^r L — k 2 and £ kj+ (u), £k,-(w) are undetermined coefficients. The 
negative sign in front of the second term is just a convention. 

The unique solution under specified initial/boundary conditions are obtained by re- 
lating the solutions on the two sides of the cavity. We designate the first integral in 



Eq. (A8) as the right-going field £k,+(z,t) and the second integral as the left-going field 
£k _ (z, t). We also put in labels to indicate the domain of the fields: £^,(z, t), £^_ (z, t), 
£^ + {uj), £^_{uj) are the fields and coefficients on the right side of the cavity, i.e. for 
z > 0, and ££,(z,t), ££_(z,t), ££ + (u>), £^_{u) are the fields and coefficients on the 
left side of the cavity, i.e. for z < 0. We also decompose the oscillator field as 

— e-^E cav>k (o;) (A9) 



The solution fields on the two sides of the cavity are related by noting that Eq. (A7) 
implies continuity of £ k (z,t) at z = 0: 



££+(0,0 -££_(o,f) =££+(«.*) -£t-{o,t) (Aio) 
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and a finite jump of -j^£k(z, t) across z = 0: 



d_ 

dz 



[E? >+ (z,t) - ££_(z,t 



d_ 

=°+ ~ d~z 



[£i + (z,t)-££ t _(z,t)] z ^ _ 



ht c d 2 
e c 2 dt 2 



E, 



cav,k 



(t) 



(All) 



Using the expansions Eqs. (A8) and (A9), we get for each frequency component, 



.ht c u 2 ~ 

•^ Ec - k(w) 



(A12) 



We again invoke the approximation k 2 « =>■ A; 2 ~ un s /c on the right hand side of 

-zwE caVjk (u;) 



Eq. (A12), which becomes 



fife 

eocn s 



Taking the inverse Fourier transform 



back to the time domain, we get 



ht d 

s£+(o,t) + ££_(p,t) -££ >+ (o,t) -££_(o,t) = -^^E cav , k 

' ' e cn s at 



(t) 



(A13) 



Eqs. (A10) and (A 13) provide two conditions which allow to express the fields on one 



side of the cavity in terms of those on the other side. The result is 



£* (0,t) 



1 



4+(M \ nt c a m 



(A14) 



In our setting, all the fields have the same circular polarization, say +. The (+) 
polarization vectors for all the slightly obliquely directed beams are of course not exactly 
parallel to each other, but consistent with the small-fc approximations we have made so 
far, we ignore these slight misalignments and use a common unit vector e + = (x+iy) / a/2. 
Then, using the notations in Ref. El we write ££ + (0,t) = C 2 E^\ nc {t)e +1 ££_(Q,t) = 
C 2 E Krca (t)e + , 8% + (0,t) = C 2 E Ktrans (t)e + , 8%_{0,t) = 0, and E cav , k (t) = £ 2 £ k (t)e+. 
(We reiterate that i? k (t) denotes the cavity oscillator field and is not to be confused 



with £ k (t), which denotes an macroscopic field mode.) From Eq. (A14), we obtain 



-E'k.transW 
-Ek,rcfl(£) 



ht r 8 



-Mt) 



(A15) 
(A16) 



2eocn s dt 

It is also clear that E^ tr!ins (t) is the radiation field at the position of the cavity and thus is 
equal to the driving radiation field in the equation for i? k (i), Eq. (TH): E^ inc = E kttians (t). 
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In the remainder of this appendix, we show that the postulated equation of motion 
for the cavity field E^(t), Eq. ([7]), leads, in conjunction with the equations for the 
macroscopic field and the QW exciton field, naturally to the interpretation of l-E^)) 2 as 
the number of photons in mode k per unit area. Our argument is based on considering 
the energy balance in our theory. From the Maxwell equations and the constitutive 



equations Eqs. (Al) and (A2), we obtain in the usual way the continuity equation for 



the macroscopic field's energy flux: 

d 



dt 



+ V ■ S = +2ht c 5(z)Re 



d_ 

dt 



E, 



where 



U e , 



nle £* ■ £ + — B* ■ B 



2Re 



Ho 

£* x — B 
Ho 



(A17) 

(A18) 
(A19) 



u em is the averaged light field energy density, including the energy of interaction with the 
background dielectric, and S, the Poynting vector, is the averaged energy flux density. 
Integrating over an infinitesimal interval across the cavity at z = gives 



2Re 



£* x — B 
Ho J 



-2Re 



2 = 4 



£* x — B 
Ho J 



ht r 2Re 



z=0_ 



£*(z 



0) • ^E cav 



(A20) 



Since the left hand side of Eq. (A20) is the net field energy flux density coming out 



of the cavity, the right hand side should be interpreted as the negative of the rate of 
energy change per unit area inside the cavity. If we integrate the right hand side over 
the cavity's plane and specialize to our setting, we get 



ht c / dxdy2Ke 



£*(z 



0) • ^E cav 



d 



ht c C 2 2Re 



<9 

F — P 

-^k^rans ^ k 



(A21) 



We now use the equation of motion of E k (t), Eq. (J7j), to eliminate -Ek jtrans , and we get 
for the right hand side of Eq. ( |A21 ) : 



d_ 
dt 



Etf - 2Re ( Q^t^Eu 



(A22) 
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fik being real. From the first term, one can see that |-Ek| 2 is to be interpreted as 
the number of photons per unit area in the cavity mode k. The second term represents 
energy exchange with the exciton field. We can cast it into a more physically meaningful 
form by using p^s equation of motion. For simplicity we illustrate this point with the 



linearized form of Eq. (10). We first write the term in question as 

' ■ ----2Re(n kP iE k ) + 2Reltt k E k -pl 



- 2Re ( n k p* k -E k 



(A23) 



For the second term, we eliminate E^ using Eq. (10) without the nonlinear terms and 



finally obtain for Eq. (A22) 



K pi 




+ 27^ Im 



' 9 ; 



(A24) 



The first term is the rate of change of energy of the linear polariton field inside the 
cavity, and the second term is the rate of energy loss by dissipation. 



We can also reduce the left hand side of Eq. (A20) to our case where all fields are in 



the same circular polarization (spin) state. We state the result here: 



n b 



£ 2 J2 2n &0 [|£k,trans(t)| 2 + |£ kj refl(*)| 2 - |^k,inc(*) H 



(A25) 



It is of course also simple to verify the equality between Eq. (A25) and Eq. (A22) by 
directly applying Eqs. (j5])-(|8]). 
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